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Abstract 


This report describes the modification of a new theory of the reflectance 
of particulate media so as to apply it to analysis of the infrared spectra 
obtained by the IRIS instrument on Mariner 9. With the aid of this theory 
and the optical constants of muscovite mica, quartz, andesite, anorthosite, 
diopside pyroxenite, and dunite, we have made modeling calculations so as 
to refine previous estimates of the mineralogical composition of the Martian 
dust particles. These calculations suggest that a feldspar rich mixture is 
a very likely composition for the dust particles. The optical constants 
used for anorthosite and diopside. pyroxenite were derived during this pro- 
gram from reflectance measurements made by F. Blaine of NASA Goddard Space 
Flight Center. Those for the mica were derived from literature reflectance 
data. Finally, a computer program was written to invert the measured ra- 
diance data so as to obtain the absorption coefficient spectrum which should 
then be independent of the temperature profile and gaseous component effects. 
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1 . INTRODUCTION 




The infrared spectra obtained by the IRIS experiment on Mariner 9 contain 
extensive data concerning the extraordinary dust storm raging on Mars at 
the beginning of the encounter. The principal features arising from the 

dust are rather broad diffuse bands centered near 480 cm ^ and 1085 cm ^ . 

\ 

Such bands are commonly observed in silicate minerals and represent the 
Si-0 bending and stretching vibrations, respectively. Their precise 
shapes and locations provide information concerning the mineral compo- 
sition and particle size distribution of the dust cloud. However, owing 
to the rather featureless bands, it has been difficult to interpret the 
spectra of the du=t in very specific ways. Comparisons with Mie theory 
calculations for quartz'*'- suggest particle diameters between 2 and 20 ym 

and the location of the Si-0 stretching band permits the Si0 ? content of 

2 3 L 

the dust to be estimated as 60 ±10%. Hunt, et al., suggest the clay 

mineral montmorillonite as a major component of the Martian dust cloud. 

They base this conclusion on the comparison of the Mariner radiance data 

with laboratory transmission spectra in which montmorillonite shows two 

rather broad bands at approximately the correct positions and lacks other 

strong features in the mid-infrared. While other published spectra of 

4 5 6 

montmorillonite are somewhat different in these respects, * ’ we believe 
that a proper comparison requires that either radiance, transmission or 
absorption coefficients be compared for laboratory data and the Martian 
experiment to eliminate shifts and other changes caused by the physical 
conditions of the measurement. For instance, it is well known that the 
bands in reflectance spectra are often shifted to longer wavelengths 
with respect to their counterparts in transmittance spectra owing to the 
influence of the real part of the refractive index. 

As the dust particles are unlikely to be spherically shaped, the Mie 

theory may not be strictly applicable. The authors have developed a 

theory of scattering by particulate media where the particles may be 
7 8 

arbitrarily shaped. ’ During this contract, our theory has been 


1 


Arthur D Little 


4 


' s ! 


.Qa.: ! 


modified to apply to atmospheric particles rather than those packed 
together on the ground. The revised theory has then been used together 
with the optical constants of typical candidate minerals and rocks to 
model the radiance spectrum of Mars. As literature values of the op- 
tical constants are scarce for mineral species, we have applied classical 
dispersion theory to measured spectra of muscovite mica, diopside pyrox- 
enite and anorthosite so as to be able to adequately cover the possible 
range of structural types for silicate minerals. The anorthosite and 
diopside pyroxenite are, of course, rocks rather than minerals but are 
representative of essentially single mineral species. 

Finally, we have written computer programs to invert the Martian radiance 
data with the aid of the known (C0 2 band derived) atmospheric temperature 
profile so as to be able to derive an absorption coefficient for the 
Martian dust that is free of the temperature profile and the Martian gas 
absorption. Such absorption coefficient spectra are perhaps the best 
vehicles for making detailed comparisons between the Mariner data and 
laboratory spectra. 
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II. THEORETICAL METHOD 

To calculate the spectral radiance of the dust cloud on Mars for any 
proposed mixture of mineral particles one must first derive the scat- 
tering and absorption cross-sections of the individual particles, then 
compute the scattering and absorption coefficients of the particulate 
medium, and finally apply a radiative transfer method to determine the 
radiance of the cloud as a function of frequency. 

To calculate the cross-sections of a particle of a given candidate 
mineral, we use modified geometrical optics for wavelengths small 
compared with the particle size (coarse particle theory) , and modified 
Rayleigh theory for wavelengths large c spared with the particle size 
(fine particle theory). 

In the coarse particle theory we consider the particle to be a polyhedron 

but treat the edges of the facets of the polyhedron as though they were 

perturbations of the smooth surface of a sphere of the same volume as 

the actual particle. Therefore, we first calculate the scattering and 

absorption by the smooth sphere by ray tracing and then adjust the results 

to take account of the effect of the edges. We do this by regarding each 

edge as an elongated ellipsoid of particle material adhering to the sphere. 

We allow the ellipsoids to have a wide range of cross-sectional shapes and 

calculate their effect statistically. It is found that edges treated in 

this way can cause a large increase in the absorption, notably in strong 
7 8 9 

reststrahlen bands. ’ ’ Some of the concepts presented in these references 
such as the contact factor and the discrete layer model apply to packed 
mineral powders and are, of course, omitted in the present application 
to a dilute suspension of particles in an atmosphere. We also omit here 
the diffuse reflection from surface asperities on the particle and con- 
sider only the Fresnel reflection from a smooth surface. In the case of 
thin mica flakes, where the perturbed sphere is clearly not a good approxi- 
mation to the particle shape, we do the ray tracing directly and average 
over all orientations of the flake. This removes any refractive scattering 


T 
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from these calculations. It can be shown that the average geometrical 
cross-section of any convex randomized particle is the total surface area 
of the particle divided by The average slant path t through a particle 
ignoring refraction is given, therefore, by 


where V is the volume and S the surface area of the particle. For the case 
of a disk of diameter D and thickness H, (1) becomes 

2DH 

t = D + 2H * 

Allowing for refraction, the average path can be shown to be approximately 


^'n 2 t 2 - t 2 + H 2 

where n is the refractive index. In our computer program we have arbi- 
trarily chosen an aspect ratio D/H = 5/1 for mica flakes. Equation {2' 
then becomes 


ilOOn 2 - 51 


In the fine particle theory where the wavelength is large compared with 
the particle size, the edges of the particle become much less important 
than the general shape of the particle. We therefore represent an arbi- 
trary particle by a smooth ellipsoid. Since the particles may have a 
wide variety of shapes, we allow the ellipsoids to vary considerably in 
their axial ratios. For our purpose it is most useful to specify an 
ellipsoid by its depolarization factors L, M and N, whose sum is 4nr, 
since these factors cause widely differing absorption even in particles 
of the same volume. The average power absorbed by a randomly oriented 
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ellipsoid of volume V and complex index m, placed in an electric field E 
with angular frequency u, is^ 


uvIe ! 2 


Im (m 2 


—1 + — •+ i 

+ L(m 2 - 1 ) . Aij + M(m 2 - 1) 4it + N(m 2 - . 1 ) 


where jto means "imaginary part of." 

This expression demonstrates clearly that, as stated above, the absorbed power 
is not independent of the shape of the small particle, as is often assumed. 
Independence of shape is true only when Im 2 - lj<<l, which is certainly not 
the case for a strong reststrahlen band. In such a region of the spectrum 
m can have values that approach the poles of the expression in square brackets 
in (5). For a sphere, L‘=M=N=4 tt/ 3, and there is therefore a single pole in m- 
space on the imaginary axis at the point m D -i/2. In a strong reststrahlen 
band the trajectory of m In m-space, as frequency is changed, runs almost 
parallel to the imaginary axis, and slightly displaced from it. In the case 
of a sphere, the pole at -i/2 therefore produces a sharp absorption peak. On 
the other hand, for a distribution of ellipsoidal particles with a continuous 
range of values L,M,N (subject to L+M+N = 4tt) , a continuous distribution of 
poles will be spread along the imaginary axis in m-space. The trajectory of 
m therefore remains close to poles over a considerable portion of the trajec- 
tory, thereby giving rise to a broad absorption maximum. 


On averaging Equation (5) over all allowable combinations of L,H,N, one 
9 

finds for the average absorption cross-section of small particles with a 
wide variety of shapes 

2n 2 vd 3 / m 2 ?.nm 2 \ 

o = Im { 

3 ~ \ n 2 - 1 / 


where v is the frequency in cm and Tid 3 /6 ° V. It is found that this 
expression does indeed give a. broad absorption band in a strong reststrahlen 
region, while a Mie theory calculation for spheres^" of the same volume V 
gives a sharp absorption band. 






It is worth noting that the validity of expressions (5) and (6) demands 
a more stringent requirement on the particle size d than the criteria 
nd/A<<l and |m|ird/X<<l given by van de Hulst^ for the Rayleigh region. 
The additional restriction is clearly related to the minimum separation 
of the m-trajectory from the imaginary axis, but has not been formulated 
owing to the lack of an explicit exact theory for ellipsoids. 

Reference 9 also gives an expression for the average scattering cross- 
section for the distribution of small ellipsoidal particles. 


The absorption mechanism just discussed acts also at the edges of large 
particles, regarded as elongated ellipsoids, and is responsible for the 
enhanced absorption mentioned earlier with respect to the coarse-particle 
theory. 

In the intermediate range of wavelengths where X~d, we estimate the 

9 

cross-sections by means of a suitable bridging formula that gives a 
weighted average of the results from the fine and coarse particle 
theories in such a way that the two theories merge slowly. 
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The radiative transfer method by which the radiance of the dust cloud 
is calculated is based on a six-stream model, also used by Conel, ^ with 
three upward-directed streams and three downward-directed streams, all 
inclined at an angle of cos ^(l//3) to the vertical. Additional programs 
were written to carry out the calculation of a brightness temperature 
based on the radiance obtained from an arbitrary number of isothermal 
layers for a dust cloud, treating the ground as a blackbody. We assume 
that the dust cloud varies with height Z only with respect to concentration 
N (particles/cm^) and not with respect to particle size distribution or 
mixing ratios of different minerals. The details are given in Reference 9. 


III. OPTICAL CONSTANTS OF MINERALS 


In order to model the radiance Spectrum of the Martian dust, it is necessary 
to have information on the complex refractive index, m = n-i k, of the can- 
didate species over the spectral range involved. As such data has long 

12 

been available for quartz, this mineral is commonly used for model cal- 
culations despite the fact that it represents an extreme in the behavior 

12 

of silicate minerals. The literature data is particularly useful as it 
is given in the form of classical dispersion theory parameters so that the 
optical constants at any point in the infrared can be calculated when 
needed from a small set of oscillator parameters. 

13 

Recently, Pollack, et al., have published similar parameters for basalt 
and andesite rocks. The difficulty with using such data for rods is that 
in addition to the usual problem of orientation of the crystallites common 
to any microcrystalline material, rocks are generally mixtures of variable 
composition. Therefore, any such parameters would be appropriate only for 
the particular rock used in the measurement. However, Prof. Frondel of 
the Harvard geology department informs us that the composition of andesites 
are not as exceesively variable as for many other rocks. Andesite consists 
principally of oligoclase or a'.Jesine feldspar but other minerals are fre- 
quently present. 

There are however, some rocks of almost single mineral composition. They 
include dunite, anorthosite and pyroxenite. We believe it is worthwhile 
to model using optical constants for such rocks as they represent struc- 
tural types and the only significant variations (beside anisotropy) are 

14 

in the cation ratios. It is well known that the Si-0 stretching fre- 
quency correlates with structural type. Such modeling amounts to disre- 
garding the anisotropy involved, on the basis that it would always be 
approximately the 6ame for natural unoriented samples. A further reason 
for the use of essentially monominerallic rocks as opposed to single 
crystals of minerals is that it is often hard to obtain suitable mineral 
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samples for analysis and that, in general, three 6ets of optical ccr.stants 
would be required for each single crystal corresponding to the three prin- 
cipal mineral orientations. 

Optical constants for dunite (olivine) calculated by R. Vincent were 
obtained from R. Curran of NASA Goddard Space Flight Center. A sample 
of anorthosite from Essex County, New York, was obtained from Ward's 
Natural Science Establishment. Prof. Frondel of Harvard University 
examined our sample and noted that some crystallites were rather large 
so that averaging might be necessary. Anorthosites consist largely of 
feldspar anorthite which is the calcium-rich end member of the plagio- 
clase feldspar minerals. Our polycrystalline sample was polished for 
reflectance measurements on three sides. The plan was that any orien- 
tation effects caused by large non-random crystals within our sample 
could be examined and might be averaged out. In addition, a microcrystal- 
line diopside pyroxenite rock, whicn is essentially pure diopside, was 
obtained from Prof. Frondel and polished for a similar analysis. 


The reflectance measurements were made by F. Blaine of NASA Goddard Space 
Flight Center. A wire grid infrared polarizer was used to ascertain 
whether significant polarization effects exist. After examining the 
measurements made both with and without the polarizer for the pyroxenite 
and anorthosite, we chose suitable spectra for reduction to the optical 
constants.. 
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The optical constants of muscovite mica have been obtained from the 
published data of Vedder.^ As the author no longer has the concomitant 
data on the refractive index, we used his recorded measurements, together 
with Blaine's measurements for the anorthosite and pyroxenite rocks to 
produce classical oscillator parameters for use in our modeling calcula- 
tions. For the two rock samples a set of averaged optical constants was 
produced while for mica three sets of optical constants for the E|ja, E | | b , 




and E_[_a,b orientations -./ere obtained. The E_|_a,b Is very close to an Ejjc 
orientation so that the three possible sets of optical constants for this 
mineral are now available. 

The computer program used to generate the optical constants from, re- 
flectance data in terms of classical oscillator parameters (Lorentz 
lines) is similar to one previously used for garnet and glass. ^ The 
old program has been modified to greatly improve the convergence rate, 
and a statistical section has been added at the end to calculate the 
standard deviations of the Lorentz line parameters. A detailed dis- 
cussion of the method of regression is given in the Appendix. 

Tables 1-5 show the Lorentz line parameters and their estimated standard 
deviations, o, obtained by fitting the measured reflection spectra. The 
formula for the dielectric constant in terms of the tabulated parameters 


€ - Cq + 


S I + fy, iL. 

J v. I v. 


where S. is the line strength, the damping and the frequency for 
each resonance, j. ^ is the high frequency dielectric constant. 


Figures 1-5 show the measured spectral data (plotted as points) and 
the theoretical spectrum for comparison. It is of course possible to 
improve the fits still further by the addition of additional resonances, 
but we believe that the experimental data does not warrant additional 
effort. We recognize that a truly rigorous set of averaged optical con- 
stants for rock species cannot be obtained in this way in view of orien- 
tation-polarization effects, the imperfect nature of the surfaces and 
tue approximation involved in averaging reflectances rather than cross- 
sections. Nonetheless, we believe that the 6et of pseudo-optical con- 
stants obtained here for monominer all ic rocks are entirely adequate to 


















characterize the structural type represented by the spectral data for 
the Martian dust. 
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0.0034 
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0.0412 

0.0019 
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: • 
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0.1047 

0.0224 
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0.084 

V 

10 

348.6 

1.3 

0.0440 

0.0144 

0.137 

0.039 

• 

e o 

2.713 

o(e Q ) = 0.057 







TABLE 1. Lorentz line parameters for Mica, E J { a . 

Range of validity is 320 5- v !? 1500 cm 
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Lorentz line parameters for Mica, E||b. 
Range of validity is 320 5! v < 1500 cm • 
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o(Sj) 

.012 
.029 
.018 
.006 
.009 
.053 
.060 
.010 

e Q ■» 2.243 a(e Q ) - 0.038 

TABLE 3. Lorentz line parameters for Mica, E_[a and b. 
Range of validity is 300 S v $ 1500 cm 
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4. Lorentz line parameters for 
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anorthosite 
< 1400 cm -1 
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Lorentz line parameters for pyroxenite. 
Range of validity is 230-5 v< 1700 cm 
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IV. MODELING THE MARTIAN DUST 
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We have used the optical constants for the various minerals and rocks 

together with our theory of particulate scattering to model the 

observed infrared spectra obtained during the Martian dust storm. 

For this purpose an averaged spectrum was obtained from Dr. Conrath, 

together with the relevant temperature profile. This profile consists 

of 19 atmospheric layer heights and associated temperatures, a scale 

height derived from the 15 u carbon dioxide band, and a blackbody 

temperature for the ground. We used a linear interpolation of the 

temperature data to get an average temperature for each layer so that 

we actually used 18 atmospheric layers in our program. Some of the 

modeling calculations had previously been run with an 8 atmospheric 
2 

layer profile taUen from the revolution 20 and the spectral shape 
differences between the two profiles are minor. 

A set of computed spectra for quartz as a function of particle size 
is shown in Figure 6 together with the Martian dust. This set illus-r 
trates the trend to feature broadening as the particle size is increased, 
which has also been observed in Mie scattering calculations.^ Thin trend 
appears to be quite general and indicates the approximate particle size 
by the breadth of the Si-0 stretching band. We treat quartz as if it 
were randomly oriented so that we always take a 2:1 ratio for the 
eJ_c:E||c particle densities. For submicron particles, our calculations, 
which are in this case based on Rayleigh scattering by a distribution of 
elliptical particles, clearly indicate that the sharpening no longer 
occurs. The absorption coefficient decreases with Nd^ and the scattering 
coefficient with Nd 6 , but by the time one reaches submicron sizes, the 
scattering is negligible compared to the absorption. Therefore, on the 
assumption that the volume fraction ud^N/6 of the particles is constant, . 
the absorption coefficient is also constant and one can no longer estimate 
N and d independently from a given spectral band. We have demonstrated 
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this difficulty in several runs and believe the Mie theory would show 
the same results for a distribution of shapes, as it passes over to 
Rayleigh scattering for very small particles, if the Mie theory could 
be extended co such a case. Owing to the calculated breadth of our fire 
particle bands, other constraints must be applied to set a lower boundary 
on particle size if the particle size is quite small. For coarse par- 
ticles, however, the broadening trend is so strong that the feldspars, 
which have smaller frequency dependent excursions in the optical con- 
stants, produce relatively featureless spectra at particle sizes in 
excess of 10 p. 

Figure 7 shows the spectra of 1 p particles of the species that we have 
run (in the case of mica flakes, 4 p diameter particles having a thick- 
ness of 0.8 p) . The general trend of the Si-0 stretching bands near 
10 p with structure can be seen to run from quart; at the high fre- 
quency end through the feldspars (anorthosite and andesite), mica, and 
pyroxenite to dunite. This trend, as is well known, is due to the 
variation in silicate structure where quartz represents a complete sharing 
of the tetrahedral SiO^ oxygens and the sharing is gradually reduced until 
in the olivine structure possessed by the dunite there is no oxygen sharing. 

The SiO, tetrahedra are, in this case, isolated by the intervening cations 
** ++ 

(mostly Mg ). The trend is of a general sort, however, and individual 
variations within each structural type are known to cause a range of band 

positions that allows considerable overlap between adjacent structural 

„ 14 

types . 

The particle densities shown in Figure 7 were held constant. They are 
adjustable parameters in our computer program but their basic effect is 
simply to change the spectral level provided the densities are sufficient 
that the spectral bands may be clearly observed. The mica data was cal- 
culated for a random orientation of flakes by taking equal numbers of 
particles of each of the three orientations Ej|a, E| |b, and E_jja,b (essen- 
tially^ E||c). We have also examined the effect of considering only 
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the two orientations.. E||a and F. J ] b » as opposed to all three orientations 
in order to examine the range of differences to be expected if the par- 
ticles are settling slowly with the sheets aligned parallel to the ground 
ather than tumbling in random fashion, although even in such a case sig- 
nificant amounts of the E||c orientation would be seen owing to the diffuse 
nature of the radiative transfer. 

We chose to obtain the optical constants and run the mica data, not only to fill 

in the structural sequence discussed above, but as a proxy for the clay 

3 

mineral mantmorillonite suggested by Hunt , et al . , to be a good Martian 
candidate. Transmission spectra published by Moenke' 5 indicate that mont- 
morillonite and muscovite have very similar spectra, which is not unex- 
pected owing to their common sheet-like structure. Mocnke shows a Si-0 
stretching frequency about 10cm ^ lower in muscovite than in montmorillonite , 
although the spectra of both species are known to be somewhat variable.^ 

As already stated, andesites are rocks of somewhat variable composition. 

The evidence that this high concentration feldspar containing rock had 
similarities to the Martian data led us originally to propose the all 
feldspar rock anorthosite as a candidate species and so to develop its 
optical constants. The resulting spectrum gave, for the first time, to 
our knowledge a possible explanation of the hitherto mysterious band near 
1380 cm However, a reexamination of the data (Figure 4) which resulted 
in the resonance giving rise to this feature casts some doubt on its 
reality. The mathematical analysis (Table 4) appears to validate the 
line in question, but the raw data shows the reflectance level to flatten 
out at still higher frequencies. Because of the low value of the reflec- 
tance and the noise characteristics of the measurements, we are hesitant 
to claim the reality of the 1380 cm feature. 

Examination of the data in Figure 7 indicates that the Si-0 stretching 
band of the two feldspar containing rocks and the muscovite mica fit the 
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Martian data considerably better than the other species with respect to 
both band shape and position. The feldspars appear to have a somewhat 
higher stretching frequency band than the mica, as expected, and thus 
fit the Martian data somewhat better. However, the low frequency band 
near 500 cm ^ is another matter. For this band the mica is preferable 
owing to its lack of a window region in the center of the band, but 
contraindicated by the relative strength of the oand compared to the 
Si-0 stretching band. In the latter respect, the feldspars are pre- 
ferable. Both of these characteristics might be Improved on by the 
choice of a suitable mixture, as the attempt to model the spectrum of 
Mars with a single species is likely to be quite unrealistic. 

The choice of 1 p for the anorthosite size in Figure 7 was based on the 
band breadth behavior of a set of anorthosite spectra analogous to 
the quartz set given in Figure 6. The Si-0 stretching band, as shown 
in Figure 7, is slightly less broad than in the Martian data suggesting 
that the particle size is slightly greater than 1 p. Our 2 p anorthosite 
results on the other hand give a broader band whose shape has been suf- 
ficiently distorted so that it no longer resembles the Martian band. 
Similarly, for mica (whose band position makes it also a tenable can- 
didate) the best appearing spectrum results from a flake of thickness 
near 1 p and significantly smaller particles have a distinctly narrower 
band. 

If we restrict ourselves to single species, we feel that in view of the 
lack of optical constants for very many mineral species, a reasonably 
good fit to the Martian data is shown by both anorthosite and mica in 
Figure 8. Here we have simply reduced the particle densities so as to 
better approximate the data. While the anorthosite fit Is by no means 
perfect, the band positions for both major bands (ignoring the 500 cm ^ 
peak in the middle of the low frequency band) , their breadths and rela- 
tive intensities strongly suggest a feldspar as a principal candidate 
for the. Martian du3t. The mica fit Is not quite as good as discussed 
above . 
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All the plagioclase feldspars appear to show similar spectral structure^ 
to the anorthosite and, as the window in the middle of the 500 cm * band is 
the principal difficulty, we decided to try a mixture of anorthosite with 
quartz. The result of one run is shown in Figure 8. The quartz does 
fill the window in the 500 cm ^ band as well as shift the Si-0 stretching 
band to slightly higher frequencies which is desirable, but, unfortunately, 
introduces a discontinuity in the Si-0 stretching band. This feature 
is due to the use of 1 p quartz rather than a slightly larger size (see 
Figure 6) to avoid a defect in our simple bridging relationship which was 
designed for mixing either coarse particles of several types or fine par- 
ticles of several types. The defect results from bridging between the two 
theories for the entire mixture rather than separately for each mineral 
component. Judging by our various runs, the best fit should be produced 
by fine particle anorthosite and somewhat coarser particle quartz. We have 
dashed in the expected result of such reprogramming on the figure. As the 
evidence of our various runs suggested that a combination of coarse particle 
(4 p) mica and quartz might also result in a better fit than the mica alone, 
we ran such a mix as well, and it is also included in Figure 8. It is im- 
portant to note that 4 p mica, as before, refers to the disk diameter and 
the thickness is taken to be 0.8 p. The relative coarseness (4 p) of the 
quartz particles is used principally to assure that the mica spectrum is 
not forced into the fine particle regime owing to the previously mentioned 
defect in our simple bridging relationship. The Si-0 stretching band is 
better centered than for the mica alone but appears too broad. Significant 
additional structure in the 800 cm 1 region and in the 500 cm 1 band is 
created but all things considered the fit is not too bad. 

In comparison to the mixture of anorthosite and quartz, the relative flat 
ness and extra breadth of the Si-0 stretching band and the intensity 
reversal shown compared to the two principal Martian bands (near 500 cm ^ 
and 1000 cm are the most significant reasons for our preference of the 
feldspar dominated mix. This is even discounting the somewhat uncertain 
1380 cm ^ feature as discussed previously. 
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We believe that a feldspar dominated mixture containing some quartz 
and mica would provide a still better fit, but this would require some 
reprogramming of the relationship bridging our fine and coarse particle 
theories. 


cu 
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V. DETERMINATION OF K(v), S(v) AND SCALE HEIGHT H OF THE DUST CLOUD BY 
INVERSION OF OBSERVED SPECTRA 

A. INTRODUCTION 

AC our meedng at Goddard Space Flight Center in July 1973 we suggested 
the posslblity that the effect of atmospheric gas absorption could be 
removed from the dust cloud spectra if, with the help of a simple radiative 
transfer theory, one suitably combined the spectra observed for two quite 
different temperature profiles. We have pursued this idea further and 
now believe that the two-temperature profile scheme, besides eliminating 
the effect of gas absorption, can provide spectral data directly related 
to the scattering and absorption cross-sections of the dust particles 
themselves and therefore be of greater diagnostic value than the cloud 
spectral radiance or brightness temperature. The parameters that can be 
derived by the two-temperature profile method are the scattering and 
absorption coefficients S(v) and K(v) of the dust cloud which are pro- 
portional to the scattering and absorption cross-sections of the particles, 
and independent of temperature distribution and gas absorption. The method 
requires an assumption about the scale height of the dust cloud, e.g., the 
same as that of the atmosphere. 

In the spectral regions with negligible gas absorption the method is very 
simple. With gas absorption a more involved computer program is necessary. 
In both cases the procedure can be readily used in computer moceling, i.e., 
in calculating spectral radiance for given temperature profile, particle 
density distribution, and particle scattering and absorption cross-sections 
derived from our theory or the Mie theory. 

The foregoing inversion scheme was worked out at a time when it was thought 
that the dust particles were large enough, relative to the wavelength. 
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that S and K would be of comparable magnitude. It has since developed, 
from aerodynamic considerations and the spectral evidence discussed above, 
that the particle size (M. u) is small compared with the wavelength 
('vlO w) « Under these conditions the effect of S on the spectral shape 
Is too small to permit S to be determined by the two-profile method. We 
have therefore considered the possibility of determining the scale height 
H and K(v) for the dust cloud, instead of S(v) and K(v), by the two- 
profile scheme. 

B. DETERMINATION OF S AND K IN SPECTRAL REGIONS HAVING NEGLIGIBLE GAS 
ABSORPTION 

We assume that the dust cloud varies with height Z only with respect to 

3 

concentration (particles/cm ) and not with respect to particle size dis- 
tribution or mixing ratios of different minerals. Under these conditions, 

S and K are each proportional to concentration n(Z). The radiance of the 
cloud, in spectral regions with no gas absorption, is not changed if the 
cloud is compressed to uniform particle concentration, provided that each 
plane of particles retains its original temperature. T he new height Z of 
any plane is related to the original height Z Q by the formula 

1 /* Z ° 

2 “ n7 / n( V dZ 0 < 8 > 

° 0 

where n^ is the concentration near the ground. 

We represent the radiation distribution at any point in the uniform cloud 
by six mutually perpendicular discrete beams equally inclined to the Z-axis. 
From symmetry the three upward-directed beams are of equal intensity and 
produce a resultant flux density I(v). The three downward-directed beams 
have a resultant flux density J(v). The radiative transfer equations then 
have the simple form 
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TV = -(K + S)I + SJ + KB 

04 

- ~ “ -(K + S)J + SI + KB 
at 

where B ■» B(T,v) is the Planck function at any point. 


(9) 

( 10 ) 


Using a Green's function approach, we have solved (9) and (10) for the 
radiative flux N leaving the top of the cloud under the assumption that 
the ground is composed of the same particles as the cloud and has a uniform 
temperature corresponding to a Planck function B Q (v). 

The result is 


N <* (1 - R)y 



H 

B(h)e“ Yh dh + (1 - R)B 0 e _Yh 


(ID 


where H is the thickness of the uniform cloud and h = H-Z. The parameters 
R and y are the reflectance (albedo) of a semi- inf inite cloud, and the 
extinction coefficient of the radiation in the cloud, respectively, and 
are related to K and S by the formulas 



y = 



■r 2KS 


( 12 ) 

(13) 


In practice one could divide the cloud into a number of isothermal layers 
as in Figure 9. Equation (11) then becomes 
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N - (1 - R)|b 4 + (B 3 - B 4 )e' Yd 4 + (B 2 - B 3 )e“ Y(d 4 + d 3 ) 

+ (B l - B 2 )e~ Y<d 4 + d 3 + + (Bq - B 1 )e" Y(d 4 + d 3 + d 2 + V} (14) 

This simple result can obviously be extended to any number of isothermal 
layers. 

Eq (14) can be modified to replace layer thicknesses by layer heights and 
to treat the ground as a blackbody rather than as an isothermal semi- 
infinite extension of the cloud. The new equation, generalized to n layers, 
is 






Qo 

<9,. v>; 


where yH “ 0, B„ is the Planck function for the ground and the H. 

O 1 

represents the heights of the i layers. The calculations have been made 
on the assumption of an exponential dependence of particle concentration 
on height. As before the cloud is treated as if it were compressed to 
layers of equal particle density having, however, the temperatures of the 
original layers. 

To determine R and y from the Martian spectra, one would select two 
spectra with widely different temperature profiles taken closely enough 
in time and space so that the dust cloud can be assumed to be unchanged. 

For the two spectra Equation (15) can be written in the form 

(16) 
(17) 


N = (1 - R) F(B a , B 3 , B 2 , B 1 , B q , y) 

N‘ = (1 - R) F(B' 4 , B’ 3 , B ' 2 , B' r B' 0 , y) 


where N. N', and the B's and B' 's are known. 


On dividing (16) by (17), 1 - R cancels out and we obtain an equation fo" 
Y alone. This can be solved readily by Newton's method. R can then be 
determined directly from Equation (16). 


The scattering and absorption coefficients can be found by means of the 
relations 


_2lR_ 

(18) 

1 - R 2 


y(1 - R) 
1 + R 

(19) 


which one can derive from (12) and (13). 
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C. DETERMINATION OF S AND K IN THE PRESENCE OF GAS ABSORPTION 


In spectral regions where gas absorption is significant, K/S is, in 
general, no longer independent of Z since the absorption coefficients of 
the gas and dust do not vary in the same way with Z. The simple solution 
(11) therefore does not apply. We may, however, make the dust cloud 
uniform as before by the Z-transformation (8) . At the same time we must 
similarly adjust the gas absorption coefficient from K (Z^) to K (Z) by 
the formula 


K g (Z) 


n 0 

n(Z Q ) 


W 


( 20 ) 


Next we divide the cloud into a number of isothermal layers as before. 
In any given isothermal layer, the general solution of (9) and (10) is 


I = Ae" YZ + Ce yZ + B (21) 

J = RAe" YZ + |ce YZ + B (22) 


where A and C are arbitrary constants to be determined by the boundary 
conditions of continuity of I and J at each interface. At the top of 
the cloud J «* 0 and we assume now that the ground has an arbitrary 
emissivity e(v). For the case of the four-layer subdivision shown in 
Figure 10, we have eight arbitrary constants, A^, C^, Aj, C^, A^, C^, A^, 
and eight equations derived from the boundary conditions. These 
simultaneous equations are displayed in Table 6. It is clear how the 
matrix can be extended to the case of more layers. 
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TABLE 6 

Simultaneous Equations for 4 Layers 
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The outgoing flux density at the top of the cloud is obtained from the 
condition of continuity of I at this interface, which is 

N = A.e ^4^4 + C.e^4^4 + B. 

4 4 4 
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For a given experimental situation B(Z) and K^(Z) from Equation (20) are 

known. For any choice of and S for the cloud of particles one can 

therefore obtain R(Z) and y(Z) from S and K = + K^. Then, with a 

suitable assumed value of e, all coefficients in Table 6 are known so that 

A, and C, can be calculated by Gaussian elimination. A theoretical value 

of N can therefore be found from (23) for the particular choice of and 

S. The same can be done for an experimental run with a widely different 

temperature profile, yielding a theoretical value N'. One then manipulates 

K and S until the theoretical values of N and N' agree with the measured 
P 

values. Instead of the two-dimensional Newton's method, which is often 
troublesome, we prefer an iterative method with over-relaxation to carry 
out the computation. 


a 

v f P 

h 
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The question of how large a difference in temperature profile is required 
to give good values of K^(v) and S(v), free of the effects of temperature 
gradient and gas absorption, is best answered by actual trial on the 
computer. For this purpose one could plot a versus S curve correspond- 
ing to matching of theoretical and experimental values of N for each 
temperature profile. Each of these curves can be obtained by Newton's 
method in one dimension. The intersection of the curves determines the 

required values of K and S. The accuracy of the procedure can be eval- 
P 

uated from the crossing angle of the two curves. 


D. DETERMINATION OF K(v) AND H WHEN S(v) IS SMALL 


We now restate the problem as being to find both the absorption coefficient 

K (v) and the scale-height H of the particulate layer given two measured 
P 

radiance spectra for widely different, but known, temperature profiles 
and a known profile of the gas absorption coefficient K (v) . Again we have 

o 

two equations and two unknowus. The equations can be derived as follows. 

In the absence of scattering the upward and downward radiative fluxes, I 
and J are given by the differential equations: 

4=- = ~KI + KB C 

a L 

- 44 - -KJ + KB C 

aZ 

where B is the Planck function and K = K + K is the total absorption 

g P 

coefficient at any height Z. Since these equations are not coupled ion 
the assumption that the ground is black) and we are interested only in 
the upwards flux, Equation (25) may be disregarded. For a set of n layers, 
each with uniform K and B, the solution of (24) is: 
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I 2 = A 2 e' K 2 Z + B 2 


r if 


. -K Z , „ 

I = A e n + B 
n n n 


At Z «* 0, 1^ ■ Bq where Bg is the Planck function for the ground temperature. 
At Z ■» Z^ the boundary between the first and second layer, I = I 2> etc. 
Therefore: 


B o - B 1 


l n 


-K 1 Z 1 -K 2 K 1 
e L \ - e A 2 = B 2 - B 1 


~ K 2 Z 2 ~ K 3 Z 2 
e a 2 - e A 3 “ B 3 “ B 2 


- K 3 Z 3 " K 4 Z 3 

e A 3 “ 6 \ * B 4 - b 3 


-K ,Z , -K Z . 

n-1 n-1, „ n n-1. „ 

e A i - e A = - B 

n-1 n n 


The radiance from the top of the cloud is given by: 


N - A e _KnZn + B 
n n 


For given values of the Ks, Zs and Bs, the pet of Equations (27) can 
readily be solved sequentially for the unknown A's, starting at the first 
equation. In this way, A^ can be computed, and therefore N from Equation 
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la order to avoid overflow problems in the computer program, we replace 
A^, A^, • • • • by X^i 5^2 , • * • • , where 


V 

“1 


- k 4V 

e A, 



(29) 


Then (27) and (28) become 


-K.Z, 

x i “ e <V" V 


-(K - K )Z K,Z - K.Z, 

X 2 - e 2 1 \ - e 2 1 . “(B 2 - B 1 ) 


(30) 


and 

N = X + B (31) 

n n 


As a preliminary step in the computer program we have chosen to compress 
the atmosphere in such a way that the exponentially distributed dust with 
assumed scale height H becomes a layer of uniform density extending from 
the ground to a height H. The required transformation is: 


Z’ 



(32) 


where Z is transformed into Z 1 . The same transformation is, of course, 
applied to the atmospheric gas. The experimentally determined temperatures 
T(Z) are kept constant in the transformation, i.e.: 

T(Z’) = T(Z) (33) 
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The known absorption coefficient K (Z) of the gas is transformed according 

B 

to the relation 


K (Z) 


V z,) " dZVdZ “ K g < Z > 6 
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In the computer program, for the first temperature profile at a given 
frequency, we select a value of the dust scale height H and then by 
Newton's method find what value of the dust absorption coefficient K p 
gives a theoretical radiance in agreement with the experimental value at. 
the chosen frequency. This is repeated for a series of values of H, so 
that a graph of K p versus H can be drawn from the computer output. The 
sams is done for the other temperature profile. The intersection of the 
two Kp versus H graphs gives the solution for dust scale height and the 
absor, tion coefficient at the chosen frequency. On repeating the procedure 
at another frequency a different value of K p will, in general, be found 
but H should, of course, be the same if the experimental data are noise- 
free and free of systematic error. 


On running the program for morning and afternoon temperature profiles at 
a frequency of 1000 cm , for which the atmospheric absorption is very 
weak, we found that the two computed K p versus H curves crossed each 
other at an angle of about 90° (Figure 11) and yielded the solution 
H '*» 7 km and K p = 1.1 x 10 ^ cm The value of H is in fair agreement 
with the scale height of about 10 km for the Martian atmosphere. The 
optical thickness K p H of the dust cloud at this frequency is 0.77. 

Unfortunately, on repeating the computations at 506 cm ^ (Figure 12), we 
found that the K p versus H curve for the morning profile stayed above 
the curve for the afternoon profile for all values of H. Therefore, no 
solution could be obtained. The reason for this difficulty is that the 
morning temperature profile is too close to isothermal to provide much 
information about the duet cloud. If the profile were strictly isothermal, 
at the same temperature as the ground, the observed radiance would have a 
blackbody spectrum corresponding to the ground temperature and would 

therefore contain no information about K (v) or H. 
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To investigate the question further, we have examined the behavior of N 
versus K for the various cases. When the temperature varies monotonically 
with height (including the ground, temperature) N versus is a steep 
curve that intersects the experimental N at a well defined value of K^. 

On the other hand, when the Planck function for the ground is between the 
Planck functions for the various layers (morning profile), the curve is 
quite flat and Kp is very susceptible to small errors in the experimental 
radiance. 


We believe that the procedure for determining K and H would work very well 
if the second profile were monotonically Increasing with height or mono- 
tonically decreasing at a significantly lower rate than for the afternoon 
profile. We recommend that NASA should search for such a pair of profiles 
and carry out the necessary computations with the help of our program. 

We also recommend, if a suitable second profile is not available, that 
the data for the afternoon profile should be used alone, along with an 
assumed value of dust scale height, e.g., 10 km, to determine a complete 
K spectrum. Such a spectrum will be free of the effects of the tempera- 
ture distribution and of atmospheric absorption and will be equivalent to 
the transmission spectrum (ulthout reflection losses) obtained by a labo- 
ratory spectrometer for a single crystal of the dust material. Such a 
presentation of the Martian infrared data would allow simple averaging of 
K over many spectra for purposes of comparison with laboratory trans- 
mission spectra of various candidate minerals. 


E. COMPUTER PROGRAMS 


A listing of MARK,, the computer program we wrote to carry out the Inversion 
when S = 0 is given in Figure 13. This program was written in double 
precision fortran for the General Electric Mark III Time Sharing Service. 
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mRK III 03EDT 09'87/74 

IOOC UP TO 50 LAYERS 

110 IMPLICIT DOUBLE PRECISION <A-H> 

180 IMPLICIT DOUBLE PRECISION <0-E> 

130 DIMENSION T<50>.P<50>.ALPf 50>.Z<50>.TEHPt 5>»AA<50> 

UO i»CAPK<50> *tS0<50> 

US REAL TEMP. DX. PIP. ZR 

150 BPLCTEI «1 . I909D-1 8* VN*VN*VN/(DEXPC | .43800* *N/TE> • t -D0> 

UO ITER*! 

1 70 | PR I NT •'•READ IN VN. AVCBAD.N.H.E.TG. 4LPM.P5.TS** 

180 READC-NUM3~*999>WN. AVGRAD.N M » EM. T C. ALPM . C PC l > . T C 1 > . 1 • | , N > 

185 999 FORMATfV> 

190 PRINT.VN. AVCR AD * N » H » EM.TG.ALPH. CPC 1 ) « TC 1 ) * 1 • I . N> 

800 CAPI-l.D-5 

205 ALPH* 1 0 «D0* • ALPH 

210C TRANSFORM P TO Z 

220 C PO IS AT THE GROUND 

230 RGASM .8892D6 

232 PC0NV-9.B693D-A 

23A DO 10 I ■ I . N 

836 10 PC I 1 *P< 1 )*PCQNV 

837 Nt«NM 

83P DO 11 1 •} . Nl 

240 I 1 ALP C I > -ALPM* 15. DO* < PCI >*PC I ♦ 1 > > / < 2 *D0*D$CRTC C T« X >*TC I *1 > >/2*D0>> 

241 CHARS «3 7 ft. DO 

260 DO 60 l«i#N| 

261 1 E*0 

862 DO 70 INM*5 

263 70 TEMPCINl.O. 

264 DDX-O.DO 

865 DPIP*TC I >/PCI > 

266 65 CONTINUE 
26T OX-ODX 

268 PIP-DPIP 

269 ER»CLC1NT< l#DX.PIP.TEHP> 

270 EC I >«-RGAS*DBLE<ZR>/GHARS 

271 1E-1EM 

272 IFME.GE.N) GO TO 51 

273 DDX»PClZ*t )*PC1E> 

274 DPIP-CTC lE>*TCIE*t >> /CpC I Z > *P< l Z* 1 »> 

275 IFC1.GE.1E> CO TO 65 

276 60 CONTINUE 
280 51 IST-0 

290 DEK».005DO*CAP1 

300 CAPKD-CAP1 -DEK 

310 25 CONTINUE 

340 DO 2 1M.N1 

350 I F< I TER »CT • 1 > CO TO 53 

360 I FC l ST «GE • I 7 CO TO 53 

380 ESOCI >«H*C 1 .DO-CEXPC -Ef I >/H> > 

395 53 CAPKC-I .732051D0*ALPC1 >*DEXPCEU >/H> 

400 55 CAPK<I >-CAPKC*CAPKO 

410 2 CONTINUE 

415 CONS •CAPKC N I > • ESQ f N ! > 

420 AAC1 >»CEPL<TG> -8PLCTC | >) )*DEXPC *CQNS> 

430 DO 90 1»2*N1 

440 AACl >»AA<1 -I >*DEXPfESQCl -1 > • C CAPK Cl > -CAPK <1 - I > > > -DEXP 
450 4CZSQCI-I > *C APKC l > -CONS >• C BPL< T C 1 > > -BPLCTU -| >>> 

460 90 CONTINUE 

470 RAD«AA<N| >«8PLCT<N1>> 

475 PR I NT »R AD 

480 1F<IST.E0.1> GO TO 22 

490 DEL-RAD-AVCRAD 

500 CAPKD»C4PKD*2.D0*CEK 

510 1 ST»1 ST* I 

580 CO TO 25 

530 22 DELK»RAD- AVCRAO 

540 DFCK-CDELK-DEL>/ce.DO*DEK> 

550 CAP I «CAP ) «CDEL*DELK>/C2.D0*DFCK> 

570 PR 1 NT * I TER .CAP l 

580 ITER-ITER*! 

590 CO TO 51 
#00 END 
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It can be rather simply converted to any other fortran service by suitable 
changes in the I/O statements and replacement or deletion of the integra- 
tion subroutine CLCINT. This subroutine is used to calculate the equivalent 
heights Z(i) for the various layers, i, from the pressure and temperature 
data if the number of layers used is less than the total Lumber given. • 

It is a modified Simpson's rule integration. The calling sequence is 
ZR = CLCINT (1, DX, PIP, TEMP). 1 is an indicator for the modified Simpson's 
rule, DX is set equal to 0 for the initial call and to the desired incre- 
ment thereafter, PIP is the integrand and TDCP is a storage array. 
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Appendix : 

General Description of Fitting Process 

The first step we take is to read the graphical output of the 
spectroneter at intervals of 5-10 cm ^ and enter the data in digital 
form into a time-sharing system. After performing a linear transforma- 
tion to correct for the zero and one-hundred percent reflectance errors, 
we make a new graph of the data on a Hewlett-Packard 7200A plotter 
which is attached to the terminal used to access the time-sharing 
system. Transcription errors and doubtful data points are easily pin- 
pointed and eliminated at this juncture. We then make a final plot of 
the purified data; reproductions of this plot are used later to assist 
in the fitting process. The purified data are also listed and punched 
in cards to serve as input to the program used to carry out the least 
squares fitting. 

Our regression (least square fitting) program requires that first 
approximations to the Lorentz line parameters be supplied as input . 

We make use of the HP plotter and time-sharing system to help obtain 
the required first guesses by plotting the reflection spectrum corre- 
sponding to the guessed Lorentz lines on reproductions of the data plot. 

It is fairly easy to adjust the parameters to ameliorate a grossly erroneous 
approximation. After two or three lines have been roughly established 
in a limited band at, say, the high frequency end of the spectrum, we 
improve the fit by means of the least squares program. Then we make a 
graph of the results and guess the parameters for a few more lines to 
expand the range of frequencies covered; again the entire set is 
improved by the least squares fitting program. In this way the full 
span of the spectral data is eventually covered. We find that this 
stepwise procedure needs much less mental and computational effort than 
an attempt to guess ten to fifteen lines simultaneously to start the 
fitting process. 
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The Regression Problem and Method of Solution 

We are given the spectral reflectance of a sample as measured by a 
known instrument and we wish to find the values of a set of classical 
oscillator parameters that give the theoretical reflectance spectrum 
that best matches the measured spectrum. We solve the problem by the 
method of least squares; that is, we vary the oscillator parameters to 
minimize the root mean square difference between the measured and 
theoretical spectra. 


Let Rj be the reflectance measured at a frequency and R(v,P) be 
the theoretical reflectance of a sample with the set of oscillator 
parameters P at a frequency v. We wish to minimize the residual 

.2 


= I (R (v,.P) - R.) 
j 2 2 


(Al) 


which we can do by setting equal to zero the partial derivatives of T 
with respect to the individual parameters. Due to the complicated nature 
of the function R(v,P), however, the resulting simultaneous equations 
are highly non-linear and must be solved iteratively. 


The function R(v,P) may be broken into two parts corresponding to 
the two planes of polarization tt and s which are parallel and perpendic- 
ular, respectively, to the plane of incidence of the beam of light falling 
on the sample in the measuring instrument. We write 


R(v,P) ■» a (v) R (v,P) + a (v) R (v,P) (A 2) 

TT TT 8 8 

where a (v) and a (v) are instrumental corrections taken to be known 
n s 

functions of frequency. 


The Fresnel formulas give R (v,P) and R (v,P) in terms of the angle 

TT S 

of incidence $ and dielectric constant c = e(v,P) which characterizes 
the electrical properties of the material. Thus 


R 

TT 

R 

s 



(A3) 

(A4) 
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c cos ji - 2 
c cos 4> + Z 


(A5) 


2 - cos 
Z + cos $ 


(A6) 


Z = /e - sin^ 0 


(A7) 
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Finally, according to classical dispersion theory, the dielectric constant 
has the form of a sum of m simple resonances (Lorentz lines) : 

m S, 

(A8) 


c(v,P) 


c + I 
o 


k=l 1 + i Y. — - (— ) 2 

k \ V 


In our work, we take the set of parameters P to be made up of the following 
quantities : 


p , are seen to be 
pk 


whence, by Equation Al, we find 
I [R(v 
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(A9) 

Pik ” 

s k 

(A10) 


Ik 


P 2k 

v k 

(All) 


' 1 


P3k ° 

"7 

(A12) 

to be solved simul. aneously 

for the parameters 

3T 

0 


3p uk 

(A13) 

we find 



P) - R ] 

3R < V 1’ P > - « 

3p, t, 
vk 

(AU) 


It is evident that the quantities r^ and r g (which are amplitude ratios) 
are analytic functions of the parameters but that the intensity 

ratios R , R and R are noc, and therefore the calculation of their 

77 S 

derivatives requires some special medicine in the form of a simple 
lemma, to wit: 
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where we have written, for short, 

°k " 1+1 P 2k V " p 3k v \ < A25 > 

Equations A16 through A25 suggest clearly the few lines of Fortran coding 
needed to compute the derivatives which occur in Equation AlA. If the 
coding is systematically organized, the machine time needed to compute 
all the required derivatives is not great in comparison to the computa- 
tion of R(v,P) itself. 

The next question to be discussed is the numerical method used for 
solving the set of simultaneous equations of the type of Equation A14. 

In the past, we have used Gauss-Seidel iteration with over-relaxation, 
solving the individual equations by Newton's method. However, we were 
plagued with frequent divergences, lines whose resonance frequencies 

wandered excessively, lines coalescing with neighbors, and so on. 

19 

Others have used the method of steepest descent at first, followed by 
Newton's method as a final phase. The objection to the descent method, 
perhaps only philosophical, is that the distance traveled in the steepest 
direction must at the same time have conflicting units ascribed to it. 

To put it another way, if the units in which the problem is stated are 
changed, the sequence of steps in the descent process will suffer a 
change which is of a substantial nature. 

It is also preferable, in automatic computation, to use a unified 
procedure rather than a sequence of two or more procedures, which entails 
the decision of when to move from one phase to the next to be made, as 
well as the chore of extra coding. 

In our recent work, reported on here, we have employed a different 
technique which is a hybrid between Newton's method and descent methods. 
What we do is to calculate a certain approximation 6P to the correction 
to all the parameters P as prescribed by Newton's method but then make 
only the correction c SP, where c is a scalar selected to minimize the 
residual. In the usual form of Newton's method, a set of simultaneous 
linear equations must be solved or, equivalently, a matrix must be 



I 


inverted. When the starting approximation is poor, the matrix may be 
singular or nearly so, and difficulties in solving the equations and 
instabilities are commonly experienced. In our method, we take a matrix 
which is an approximation to the matrix occurring in Newton's method but 
which is positive definite, so that inversion is always possible. 

To elucidate the method mathematically we will eliminate the double 

subscripts p and k used hitherto to distinguish the parameters by setting 

p = p (A26) 

*£ *pk 

with £ = 0 for p = k = 0 

£ = 3(k - 1) + p for p, k £ 1 

Using this notation we may write the simultaneous equations to be 
solved (Equation A14) in the form 

1 [R(v ,P) - R ] ffiV- J L i gi « o <A27) 

i * J 3p c 


Now suppose that we have an approximation P q to the set of parameters P. 
We substitute for P in Equation A27; of course, the right-hand side 
will not be zero but some other quantity, say 6T„. Thus 


E [R(v ,P ) - R ] 3R .^ V 3l!£l = 6T. 

j J ° j 3p £ 


Now let the error in P be 6P. We substitute the expression 

o 


j i.:- 

* ! 

I 1 1 

i n 


P => P - 6P 
o 


into Equation A27 and expand in power series, keeping only the first order 
terms, to find 

ll ( 5R(v 3» P o ) 5R(v j ,Po) + [R(v ,P )-R ] 5p = 6T. (A30) 


3p 


j O j op^ 



- .1 


. ) 


Equation A30 represents a set of linear simultaneous equations to be 

solved for the errors <Sp^ and is the mathematical expression of Newton's 

method. In our procedure we drop the second tej.m ; 

2 

[R(v ,P ) - R ] .? H ^ , in the coefficient of 6p for two reasons. 

J 0 1 >, t 8p x » 

First of all, when P q is close to P, the theoretical spectrum should be 

close to the measured spectrum, so the quantity [R(v.,P ) - R ] 6p. is 

jo j A 

second order. In the second place, the matrix whose elements are 


_ 3R(V4,P 0 ) 3R(v4,P 0 ) 

a J — . — J — _ is positive definite and symmetric, and therefore 


3p 


3p 


its inverse is known to exist and is likely to be easy to compute. There- 
fore we set up and solve the equations 


T.T. • p °> aR <“.1- 1 ’°> ip . „ 

Xj 3p^ 3p^ X l 


(A31) 


Having found the elements 6p^ of the set 6P from Equation A31, we next 
find the scalar c which minimizes the residual 

T(c) ■= I[R( Vj ,P o - c6P) - Rj] 2 (A32) 


We do this simply by tabulating T(c) at a few well chosen values of c 
and finding the minimum of an interpolating parabola. Once the new 
approximation P q - c6F is found the process is repeated until the 
residual converges to some irreducible value, whereupon the last approxima- 
tion to the set of parameters is declared to be the required set. 


Statistical Analysis 

Once the Lorentz line parameters themselves are determined it is 
desirable to estimate their standard deviations to get an idea of the 
errors that may be in the parameters determined by the least squares 
fit. To do this we make the assumption that the difference between the 
theoretical and measured spectra is entirely unsystematic and due only to 
noise which is uniform over the spectrum. This may be well warranted 
in some cases but it is not entirely "alid in others, and therefore thr 
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results of the statistical analysis may not be quantitatively correct. 
The statistical estimates of error should, however, serve to Indicate 
which parameters are better or worse determined than others in the same 
spectrum. In any case, this basic assumption should always be borne 
in mind when one is interpreting the estimates of standard deviations. 

In order to calculate the standard deviations of the Lorentz line 
parameters we imagine the measurements R^ to be in error by amounts 
6Rj and calculate the change, <5P, that this causes in the parameters 
that result from the least squares fit. Evidently the equations to be 
satisfied are 

3R(va ,P + 6P) 


r [R(v ,P + 6P) - R - 6R ] 

J j i J a Po 


(A33) 


Now we expand Equation A33 in the manner used to expand Equation A28 and 
again neglect all but the lowest order terms to obtain a set of simul- 
taneous linear equations for the errors 6p^ in the parameters: 


rr 3R(v 1t P) BRCv^P) _ r 3R(vj,P) 

3p„ 3p. X j 3p c 


6R 


i 


(A3 4) 


\ 


h- 


18 v 


Equation A34 may be written in main;; notation as 
UU T 6P «* U6R 

where U is a matrix whose elements are 




aRCvj.p) 


3p i 


6R is the vector whose elements are 6R 
the superscript T denotes transpose. 

,,T 


j 


(A35) 


The matrix VJU is the symmetric, positive definite matrix encountered 
previously in the discussion of the descent method. 

Clearly, we cannot solve Equation A35 for 6P since we have no way of 
knowing 6R, but we can obtain statistical information about 6P. To start 
with we take the expectations (symbolized by the angle brackets) of both 


L ■ 
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sides of Equation A35: 

UU T <6P> ■= U <6R > (A3b) 

Obviously, if <6R> = 0, as we have a right to assume it is, then it 
follows that (6P> = 0. In other words, if the measurements are unbiased, 
then the estimates of the parameters are unbiased. To get a second 
moment we multiply both sides of Equation A35 by their transposes and 
take the expectation of the result. We may write 

UU T <6PdP T > UU T = U <6R6R T > U T (A37) 


We assume now that the measurement errors are independent and all have the 

•p 

same standard deviation a. Thus the matrix ^6R6R is a diagonal 

2 

matrix with diagonal elements all equal to a and Equation A37 becomes 

UU T <6P6P T > UU T = UU T o 2 . (A38) 

whence the matrix of variances and covariances of the Lorentz line param- 
eters is seen to be 


<6P6P T > =. (UU T ) _1 o 2 (A39) 

2 

An estimate of o may be obtained from the residual after the least squares 
fit . We take 



E[R(v ,P) - R.] 2 

J d “ 

N - 3m - 1 


(A40) 


where N is the number of measurements made 

m is the number of Lorentz lines used for fitting 


Finally it is to be noted that the variances of the classical 
oscillator parameters c , S, , y, pnd v are related to the elements of the 

O K K K. 

variance matrix computed according to Equation A39 by the relations 


< 6e o 2 > = < 6 Poo 2 > (AA1) 

<6S k 2 > = < s P lk 2 > (M2 > 
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